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ABSTRACT 


A  stiff  differential  equation  is  characterized 
by  the  fact  that  the  real  parts  of  some  of  the  eigen¬ 
values  of  the  differential  equation  system  are  negative 
and  are  very  much  larger  in  magnitude  than  the  real 
parts  of  the  other  eigenvalues. 

Classical  difference  methods  restrict  the  stepsize 
due  to  stability  conditions.  Many  authors  have  designed 
new  methods,  in  which  a  large  stepsize  can  be  taken.  We 
have  surveyed  some  of  the  important  methods  designed  for 
stiff  systems. 

Also,  we  have  proposed  an  improvement  on  the 
Liniger  and  Willoughby  method,  by  developing  a  fourth 
order  exponentially  fitted  multistep  method.  The 
computational  results  show  that  for  linear  systems  with 
constant  coefficients,  the  proposed  method's  performance 
with  a  stepsize  2h  is  the  same  as  the  performance  of  the 
Liniger  and  Willoughby  method  with  a  stepsize  h.  Even 
though  the  proposed  method  involves  more  function  evalu¬ 
ations  than  the  Liniger  and  Willoughby  method,  we 
conclude  that  it  is  more  efficient  for  linear  systems. 
For  non-linear  systems,  the  proposed  method  should  only 
be  used,  if  higher  accuracy  at  the  cost  of  more  function 
evaluations  is  desired. 
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CHAPTER  I 


INTRODUCTION 


In  physical  problems  such  as  those  occurring  in 
the  study  of  chemical  kinetics,  nuclear  reactions  and 
control  systems,  we  often  encounter  "stiff"  differential 
equations,  which  are  difficult  to  solve  even  by  numerical 
methods.  These  stiff  equations  [5]  may  be  characterized 
by  the  fact  that  the  real  parts  of  some  of  the  eigenvalues 
of  the  differential  equation  system  are  negative  and  are 
very  much  larger  in  magnitude  than  the  real  parts  of  the 
other  eigenvalues. 

Classical  difference  methods  for  solving  stiff 
systems  are  inefficient.  Stability  requirements  which 
depend  on  the  eigenvalues  of  the  system,  force  a  choice 
of  stepsize,  which  is  too  small.  Recently  however,  many 
authors  have  introduced  new  methods  [2]  designed  espe¬ 
cially  for  solving  stiff  systems.  It  is  the  primary 
objective  of  this  thesis  to  study  these  methods  and  to 
suggest  improvements  for  the  Liniger  and  Willoughby 
method  [ 19 ] . 

Chapter  I  is  intended  to  give  a  general  introduction 
to  the  problem  of  stiffness  in  ordinary  differential 
equations.  In  Section  1.1,  classical  single  step  dif¬ 
ference  methods  are  defined  and  their  shortcomings  for 
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solving  stiff  systems  are  illustrated.  In  Section  1.2, 
multistep  methods  are  defined  and  some  stability  defini¬ 
tions  to  be  used  in  later  chapters  are  given. 

Various  modified  techniques  such  as  Lawson's  Runge- 
Kutta  method  and  the  Liniger  and  Willoughby  method  are 
discussed  in  Chapter  II.  Gear's  method  and  our  proposed 
fourth  order  method  are  described  in  Chapter  III.  Finally, 
in  the  last  chapter,  results  of  test  problems  using  both 
the  Liniger  and  Willoughby  method  and  our  proposed  method 
are  given.  It  is  shown  that  for  linear  systems  with 
constant  coefficients,  the  proposed  method's  performance 
with  a  stepsize  2h  is  same  as  the  performance  of  Liniger 
and  Willoughby  third  order  method  with  a  stepsize  h.  The 
principal  advantage  of  the  proposed  method  is  that  of 
higher  accuracy.  However,  the  cost  of  function  evaluations 
is  more  than  the  Liniger  and  Willoughby  method.  The 
tradeoff  for  one  non-linear  system  is  discussed  in  the 
last  chapter. 


- 
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1 • 1  Single  Step  Methods  for  Stiff  Systems 
Definition 

Any  method  for  solving  a  differential  equation 

(1.1.1)  y  =  f  (x,y) 

in  which  the  approximation  Yn+1  to  the  solution  at  the 
point  xn+^  can  be  calculated  from  only  x^,  y^  and  h  = 
xn+1  -  xn  is  called  a  single  step  method. 

By  a  stable  method,  we  mean  one  for  which  the  error 
£n  =  y  -  y  (x  )  remains  bounded  as  n  -*■  00 .  To  see  the  pro¬ 
blem  that  arises,  when  using  single  step  methods  for 
solving  stiff  equations,  consider  the  equation 

(1.1.2)  y(x)  =  X (y  -  F  (x) )  +  F (x)  ,  A  <  0 

with  the  solution 

(1.1.3)  y(x)  =  c.exp(Xx)  +  F (x) 

where  c  is  a  parameter.  Assuming  that  F(x)  is  not  itself 
exponential,  we  see  that  the  term  c. exp (Ax)  is  significant 
for  small  x,  but  is  insignificant  for  large  x.  Neverthe¬ 
less,  the  parameter  A  appearing  in  c.exp(Ax),  because  of 
stability  considerations,  restricts  the  stepsize  both  for 
small  and  large  x. 


; 

- 


For  example,  on  applying  Euler's  method  Yn+]_  = 

Yn  +  hfn  to  (1.1.2)  with  F (x)  =  0 ,  we  obtain 

yn+l  =  V1  +  Xh)  • 

Since  y(x)  -*■  0  as  x  -*  °°,  Euler's  method  is  stable  if  and 
only  if  0 <  j  h A | < 2 .  The  behaviour  of  y^  in  the  case  of 

•  ,  .  i  i 

instability  i.e.,  j  h X | >2  is  illustrated  in  Fig.l  for 
general  F(x).  The  stability  restriction  on  h,  must  be 
maintained  for  all  x,  even  though  the  term  in  (1.1.3) 
containing  A  becomes  negligible  for  large  x.  The  ineffi¬ 
ciency  of  Euler's  method  for  solving  (1.1.2),  arises  from 
the  fact  that  for  large  x,  stability  requirements  force  a 
choice  of  stepsize,  which  is  too  small. 

This  difficulty  persists  with  higher  order  single 
step  methods,  such  as  the  Runge-Kutta  fourth  order  method. 
On  applying  the  R-K  4  method 

yn+l  =  yn+  (K0  +  2K1  +  2K2  +  K3)/6 

where 

Kq  =  hf(yn)  ,  Kx  =  hf(yn  +  Kq/2)  , 

K~  =  hf  (y  +K-./2)  and  K0  =  hf  (y  +  K0) 

2  2  n  r  3  2  n  2 

to  the  equation  y  =  Ay,  we  obtain 


t 


All  figures  are  given  at  the  end  of  the  section. 


,f  i  i  *  , 
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where 


For  stability,  we  must  have 


P (hX)  <  1 


which  gives  the  bound  0<|hX|<2.78,  when  X  is  real.  The 
stability  region  for  the  complex  X  is  shown  in  Fig.  2. 

We  now  extend  the  above  arguments  to  linear  systems 
of  differential  equations 


(1.1.4) 


y  =  Ay 


with  the  initial  conditions  y(0)  =  y  ,  where  A  is  an  n x  n 


matrix  with  n  distinct  eigenvalues.  If  the  eigenvalues 
of  A  are  in  the  negative  half  plane  (Real  X  <  0),  the  solu¬ 
tion  of  the  system  is  said  to  be  inherently  stable. 

For  h  >  0  and  x^  =  nh,  n=0,l,2,...,N,  the  solu¬ 
tion  of  (1.1.4)  is 


y(xn)  =  [exp  (hA) ] nYQ 


(1.1.5) 


Since  the  eigenvalues  of  A  are  distinct,  there  exists  a 
matrix  P  such  that 


(1.1.6) 


A  =  P  A  P 


-1 


- 

' 
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where  A  =  diag(A-^,  A2  /  •••/  A  ).  The  solution  of  (1.1.4) 
can  then  be  written  as 


(1.1.7)  y(xn)  =  P[exp(hA)]n  P  1  Y0  . 

If  either  Euler's  method  or  the  R-K  4  method  is 
used,  the  approximate  solution  is 


(1.1.8) 


yn  =  [M(hA)]n  yQ 


where 


M(hA)  = 


I  +  hA 


,2  2  ,3  3  ,4  4 

I  +  hA  +  +  ^4- 

2  6  2  4 


respectively.  With  A  given  by  (1.1.6),  equation  (1.1.8) 
becomes 


(1.1.9) 


yn  =  P[M(hA)  ]n  P  1  yQ  . 


Since  A  is  a  diagonal  matrix. 


M  (hA) 


li 


1  +  hA^ 


,22  ,3,3  ,44 
h  A.  h  A-  h  A- 

1  +  hA±  +  +  -1~- 


and  we  see  that  M(hA)^  is  a  truncated  Taylor  series 
expansion  of  exp(hA). 

Comparing  equations  (1.1.7)  and  (1.1.9),  it  is 
clear  that  Euler's  method  and  the  R-K  4  method  respectively 
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are  unstable,  if  and  only  if  for  some  i 


and 


1  +  hXi|  >  1 


h2X2  h3X3  h4X4 
1  +  hXi  +  +  +  _ 


>  1 


To  make  this  remark  clear,  we  consider  the  following 


examples 

from 

two  space: 

(1.1.10) 

• 

Y 

=  Ay 

y (0 )  =  [0.9 

(1.1.11) 

• 

z 

=  Dz 

1 — 1 

II 

0 

N 

where 

A  = 


-500.5  499.5 

499.5  -500.5 


and  D  =  diag(-l,  -1) 


The  solution  of  (1.1.10)  is 

(x)  =  exp(-x)  -  0.1  exp(-1000  x) 


Y2  (x)  =  exp(-x)  +  0.1  exp(-1000  x) 


and  that  of  (1.1.11)  is 

z1(x)  =  exp(-x) 


z  2 (x)  =  exp (-x) 
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Rewriting  y1(x)  and  y2  (x)  ,  we  have 

y1(x)  =  z1(x)  -  0.1  exp (-1000  x) 

y2 (x)  =  z2  (x)  +  0.1  exp(-1000  x)  . 

The  graph  of  the  solutions  is  given  in  Fig.  3. 

Observe  that,  except  for  those  points  in  some  neigh¬ 
bourhood  of  the  origin,  the  solutions  y(x)  and  z (x)  are 
practically  identical.  Indeed  on  most  computers,  single 
precision  computations  would  completely  ignore  the 

transient  term  exp (-1000  x)  for  x  >  0.02  (exp (-1000  x)  s 
-9 

2.06  x 10  ,  when  x  =  0.02).  We  should  therefore  expect 

that  any  particular  method,  when  applied  to  (1.1.10)  or 

(1.1.11) ,  should  be  equally  effective  at  least  for  x>  0.02. 
That  this  is  not  the  case  with  classical  difference 

methods,  is  a  consequence  of  the  great  disparity  in  the 
eigenvalues  of  the  two  systems  (A^  =  X2=  -1  in  (1.1.11), 
while  A^=  -1,  X2=  “1000  in  (1.1.10)).  Thus,  although  the 
transient  term  exp (-1000  x)  is  negligible  for  x>  0.02,  it 
is  precisely  this  term,  obtained  from  the  eigenvalue 
X2  =  -1000,  which  determines  the  maximum  permissible 
stepsize  for  a  stable  difference  scheme.  For  example,  if 
Euler's  method  is  applied  to  (1.1.11),  the  stability  con¬ 
dition  is 

(1.1.12)  max|l  +  A^h |  =  | 1  —  h |  <  1  , 


■ 

. 
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which  is  satisfied  if  |h|  <  2,  but  if  it  is  applied  to 
(1.1.10),  the  stability  condition  is 

(1.1.13)  max|l  +  X^hl  =  |l  -  1000  h|  <  1  , 

which  is  satisfied  if  |h|  <  2  xlO  .  When  solving  (1.1.10) 
by  Euler's  method,  restriction  (1.1.13)  on  the  stepsize 
may  indeed  be  reasonable,  since  in  this  interval,  the 
eigenvalue  X^=  -1000  is  certainly  significant.  Outside 
this  interval  however,  equation  (1.1.10)  is  essentially 
the  same  as  (1.1.11)  and  the  more  reasonable  restriction 
(1.1.12)  on  the  stepsize  should  be  imposed.  In  fact, 
methods  designed  for  stiff  systems,  ignore  the  transient 
terms,  as  the  computation  proceeds.  This  fact  is  what 
distinguishes  them  from  the  classical  methods. 


. 


I  1 1  '  1 
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Fig. 2  Upper  Stability  Region  for  R-K  4  Method 


Fig. 3  The  Behaviour  of  y  and  z 
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1 . 2  Multistep  Methods  and  Stability 


Define  the  multistep  method  for  (1.1.1)  by  the 
equation 


(1.2.1) 


l  (a. 

j=0 


Yn-j 


+  h  3, 


fn-i 


0 


where  a^.  and  3k  are  real  constants  with  |  |  +  |  $k  |  ^  0 

and  c*q  ^  0 .  Equation  (1.2.1)  is  commonly  expressed  as 


p(E)yn-£  h0(E)fn-k  =  0 


where 

k 

P  ( y )  =  l  a.  y 

j=0  D 


(1.2.2) 

k 

a(y)  =  l  8-  y 
j  =  0  3 


and  E  is  the  usual  shift  operator  defined  by  Eyn  =  yR+1 • 
Applying  method  (1.2.1)  to  the  single  equation 
y  =  Ay ,  yields  the  numerical  solution 


where  g^^ ,  g2  ,  .  .  .  ,  gk 
initial  values  and  the 
equation  given  by 


+ 


n  . 

g2y2  + 


n 

gkyk 


are  constants  depending  on  the 
y^.  satisfy  the  characteristic 


l  (a.  +  hX 3 ■ )  yk  j  =  0  . 
j  =  0  3  3 


(1.2.3) 


'  ■  H II  fl|§ 
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If  y ^  is  the  principal  root,  then  for  small  h 
y- ]_  =  exp(hA)  +  0  (hp+1) 

where  p  is  the  order  of  the  method.  If  X  <  0,  the  numerical 
solutions  should  tend  to  zero  as  n  00 .  Thus,  if  a  method 
is  applied  to  any  inherently  stable  system,  the  principal 
root  of  (1.2.3)  should  be  less  than  one.  Dahlquist  [6] 
has  defined  this  property  to  be  A-stability. 

Definition 

A  k-step  method  with  fixed  stepsize  h  >  0  is  A-stable, 
if  the  numerical  solution  of  any  inherently  stable  diffe¬ 
rential  equation  tends  to  zero  as  n  00 . 

If  |y^|  >  | y j |  ,  j=2,3,...,k,  the  solutions  determined 

by  the  k-1  extraneous  roots  1^'  ^3'  •••r  y^  are  small 

compared  with  the  principal  solution  g^y^.  However,  if 
| y1 1  <  | y j |  for  any  j,  then  the  method  is  unstable.  We 

note  that  for  X  <  0 ,  y(x)  -*■  0  as  x-*-  00  and  therefore  the 
roots  yk  must  be  in  the  unit  circle. 

Definition  :  Gear  [12] 

A  multistep  method  is  stable  if  the  roots  y^  of  the 
polynomial  P(y)  are  such  that  | y  ^ |  <1  and  those  roots 
which  satisfy  | y ^ |  =1  are  simple. 

Definition  :  Gear  [12] 

A  multistep  method  is  absolutely  stable  for  those 
values  of  hX ,  where  roots  of  (1.2.3)  are  <1  in  absolute 


value . 


' 
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Dahlquist  [6]  has  shown  that  the  order  p  of  an 
A-stable  multistep  method  cannot  exceed  2  and  that  the 
smallest  truncation  error  of  those  methods  of  order  2, 
is  obtained  by  the  trapezoidal  rule  (see  Fig. 4  for  the 
A-stability  domain).  Thus,  if  multistep  methods  of  order 
p  >  2  are  to  be  used,  the  constraints  imposed  by  A-stability 
must  be  relaxed. 

Several  authors  have  given  stability  definitions 
that  relax  the  requirements  of  A-stability.  These  defi¬ 
nitions  and  methods  associated  with  them  are  surveyed  in 
Bjurel  et  al .  [2] .  In  this  thesis,  in  addition  to  A- 
stability,  we  shall  use  the  concept  of  stiff  stability 
given  by  Gear  [11] . 

Definition 

A  method  is  stiffly  stable  if  it  is  absolutely 
stable  in  the  region  R-^  (Re(hA)  <  6^)  (see  Fig. 5)  and 
is  accurate  in  the  region  R^  (6-^<  Re(hA)  <  6^/ 

| Im (hA ) |  <  y ) . 

We  are  only  interested  in  stability  for  the  decay¬ 
ing  solutions  of  an  inherently  stable  system,  since  we 
choose  h  such  that  exp(hA)  is  very  small.  We  are  concerned 
both  with  accuracy  and  with  stability  for  the  other  roots. 
If  hA  =  x  +  iy  ,  Gear  shows  that  the  maximum  stable  and 
accurate  region  is  a  rectangle  bounded  by  5-^  <  x  <  S2  and 
-  Y  <  y  <  y  t  where  exp  (6^)  is  the  maximum  amount  of 


growth  allowed  in  one  step  and  y  is  chosen  so  that  the 
method  is  accurate  for  the  imaginary  components. 
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Im (hX ) 


Fig. 4  A-stability  Domain 


STABL 


Fig. 5  Stiffly  Stable  Domain 


»  Re(hX) 


^  Re (hX ) 
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CHAPTER  II 


SINGLE  STEP  METHODS 


2 . 1  Modified  Runge-Kutta  Methods 

For  the  explicit  R-K  4  method  of  section  (1.1), 
stability  requirements  for  (1.1.1)  force  a  choice  of 
stepsize  which  is  too  small.  To  increase  stability, 
several  authors  have  proposed  modifications  to  the  R-K 
methods . 

For  example,  Rosenbrock  [22]  has  formulated  the 
semi-implicit  R-K  method  for 


(2.1.1)  y  =  f (y) 


as  follows: 


(2.1.2) 


where 


K±  =  hf(yn)  +  ha1A(yn)K1  , 


and 


A (y )  =  3f (y) /3y . 
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In  his  paper,  Rosenbrock  suggests  one  possible 

choice  of  the  coefficients  a  ,  b  .  and  c  . ,  which  are 

P  PD  PD 

summarized  in  row  1  of  Table  1.  Calahan  [4]  suggests 
an  alternative  choice  of  coefficients,  which  are  summa¬ 
rized  in  rov;  2  of  Table  1.  Lapidus  et  al.  [15]  report 
that  the  method  with  Calahan' s  coefficients  is  computa¬ 
tionally  faster  than  most  of  the  single  step  methods  for 
moderately  stiff  systems.  However,  for  a  highly  stiff 
system,  the  method  is  computationally  slower  than  other 
available  methods,  such  as  the  Liniger  and  Willoughby 
method  considered  in  section  2.2.  Also,  Lawson's  method 
[16]  (see  pp.  17 )  should  be  better  for  linear  systems 
than  Rosebrock's  method  with  Calahan 's  parameters. 

Table  1 

The  Coefficients  of  Rosenbrock ' s  Method 


Method 

al  a2 

b21 

i — i 

CM 

O 

W1 

W2 

Stage 

Order 

Rosenbrock 

1-1//2 

(/2-l)/2 

0 

0 

1 

2 

2 

Calahan 

0 .78868 

-1.15470 

0 

0.75 

0.25 

2 

3 

To  provide  the  reader  with  an  intuitive  notion  of 
the  difference  between  Rosenbrock 's  method  and  the  classi¬ 
cal  R-K  methods,  we  apply  Rosenbrock ' s  method  to  y  =  Ay. 


. 


Equation  (2.1.2)  now  becomes 


y 


n+1 


R(Q)y 
2  n 


where  Q  =  hA.  Corresponding  to  Rosenbrock's  and  to 
Calahan's  choice  of  coefficients,  we  have 


R(Q) 


1  +  (/2  -  1 )  Q _ 

1  +  (/2-2)Q  +  (1.5  -  /2)Q2 

1  -  0.578Q  -  0.456Q2 

1  -  1.578Q  +  0.622Q2 


respectively . 

We  observe  therefore  that  while  classical  R-K 
methods  produce  polynomial  approximations  to  exp (Q) 

(c.f.,  section  1.1),  Rosenbrock's  method  produces 
rational  approximations  to  exp(Q).  Thus,  since  we 
might  expect  that  functions  with  transient  terms  can 
be  more  easily  approximated  by  rational  functions  than 
by  polynomials,  Rosenbrock's  method  should  be  better  for 
solving  stiff  systems  than  the  classical  R-K  methods. 

Another  modification  due  to  Lawson  of  the  classica 
R-K  method  can  be  developed  by  using  the  following  trans 
formation 


(2.1.3)  z(x)  =  exp (-xA) y (x) 

where  A  is  a  real  matrix.  The  motivation  for  using  this 
transformation  is  to  separate  a  stiff  system  into  a 


r- 
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linear  part  and  a  transient  part.  Equation  (1.1.1) 
after  transformation  becomes 


(2.1.4) 

z  (x)  =  exp (-xA)  (f (x,exp (xA) z (x) )  -  A  exp  (xA) z  (x) ) 

=  g (x, z)  . 

The  Jacobian  of  (2.1.4)  is 

3g/3z  =  exp(-xA) (3f/3y  -  A)exp(xA)  . 

Thus  the  eigenvalues  of  3g/3z  are  those  of  3f/3y-A. 
Applying  the  general  R-K  method  to  (2.1.4)  yields 


where 

P 

z  =  z  +  V  w.  K.  j=2,3,...,p 

n+1  n  .  l  i  J  ^ 

3=1 

K1  =  exp (-xnA)  (f (xn,exp (xrA) zr)  -  A  exp(xnA)zn)  , 


13 

i_j. 

II 

j-1 

= exp ( (x  +  c.h)A) (z  +  h  T  a. .K.) 

n  n  n  .  ill 

J  i=l  J 

and 

K  .  = 
3 

=  exp  (-(x  +c.h)A)(f(x  +c.h,p.)  -  Ap  .  )  . 

^  n  j  n  3  3  3 

On  transforming  to  the  y-coordinate  system,  we  obtain 


(2.1.5) 

r,  * 

yn+-^=  exp(hA)y^+  h  £  w  jexp  (  ( 1-c  ^  )  hA)  K  ^ 

where 

Ki  =  f(xn'yn>  "  Ayn  ' 

*  j“l  * 

p.  =exp(c.hA)y  +  h  l  a  .  .  exp  (  (c  . -c  .  )  hA)  K . 

J  J  n  J  1 
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and 


* 


* 


* 


for  0  =  c^  <  c2 


< _ £  c  <  1. 

P 


Lawson  proves  that  the  method  is  A-stable.  System 
(2.1.5)  is  used  computationally  for  integration  and 
Diagonal  Pade  approximants  are  used  for  approximating 


evaluated  frequently,  is  used  for  the  matrix  A.  The 
main  advantage  of  the  method  is  the  ability  of  transfor¬ 
ming  a  stiff  system  into  a  simpler  non-stiff  system, 
which  is  easy  to  integrate.  The  advantage,  however,  must 
be  weighed  against  the  disadvantage  of  evaluating  the 


Jacobian . 


' 
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2 . 2  Liniger  and  Willoughby  Method 

Liniger  and  Willoughby  [19]  developed  integration 
formulas  based  upon  the  following  identities 


(2.2.1a)  y  (x+h) -y (x)-h ( (1-y) f (x+h)+yf (x)  )  =  ex  (x)  , 

(2.2.1b)  y  (x+h)-y (x)-h( (1+a) f (x+h)  +  (1-a) f (x) ) /2 

+h2  (  (b+a) f (x+h) - (b-a) f (x) ) /4  =  e2 (x) 

and 


(2.2.1c)  y (x+h) -y (x) -h ( (1+c) f (x+h) + (1-c) f (x) ) /2 


+hz  (  (l+3c) f  (x+h)- (l-3c) f (x) )/12  E  e3 (x) 


where 


e1  (x)  =  -h4 


(0-y ) f  (x+0h) d9  , 


0 


e2(x)  =  (hJ/4) 


(20  -2  (1-a) 0  +  (b-a) )f (x+0h)d0 


0 


and 


e..  (x)  =  -  (h  /12 ) 


0 (20  —  3  (1-c) 0  +  (l-3c) )f (x+h)d0 


0 


Here  y,  a,  b  and  c  are  real  parameters.  These  formulas 
lead  naturally  to  the  following  finite  difference  schemes 
(which  we  label  and  F  ^ 2  ^  for  convenience): 


(2.2.2a) 


(1) 


y 


n+l 


=  y„  +  h ( (1-y ) f 


n 


n+l 


yf  ) 

r  n 
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(2.2.2b)  F(2)  :  yn+1=yn+  h((I+a)fn+1+  (l-a)fn)/2 

-  h2  (  (b+a)fn+1-  (b-a)fn)/4  ^ 

and 

(2.2.2c)  F(3)  :  y  , =  y  +  h((l+c)f  +  (l-c)f  )/2 

2 n+1  2 n  n+1  n  ' 

-  h2 (  ( l+3c ) f  ^  -  (l-3c) f  )/12 

n+1  n  ' 

We  note  that  formula  is  of  order  i  and  F  ^ 2 ^  is  a 

(2 ) 

special  case  of  F  with  b  =  1/3  and  c  =  a. 

In  order  to  find  a  rationale  for  choosing  the 
parameters  in  the  above  formulas,  we  follow  the  ideas 
of  Liniger  and  Willoughby  and  apply  the  formulas  to  the 
equation  y  =  Ay;  Re (A)  <  0.  The  true  solution  is  given 
by 

(2.2.3)  yn+1  =  exp(Q)yn 

where  Q  =  hA. 

Corresponding  to  equations  (2.2.2),  the  results 
obtained  in  each  case  are  of  the  form 


t2-2-4>  yn+i  = 

where 

(2.2.5a)  r (1) (Q) 

(2.2.5b)  r (2) (Q) 


(Q)y 


n 


=  (1  +  yQ) / (1  -  (l-y)Q)  , 

=  (4  +  2  (1-a)  Q  +  (b-a)  Q2 )  /  (4  -  2(l+a)Q 


. 


+  (b+a)Q2) 


. 
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and 

(2.2.5c)  r  (3)  (Q)  =  (12  +  6  (l-c)Q  +  (l-3c)Q2)/ 

(12  -  6 (1+c) Q  +  ( l+3c ) Q2 ) 

respectively. 

For  A-stability,  we  must  have  |r^^  (Q)  |  ^ 
which  gives  the  following  conditions  for 

(2.2.6a)  F(1)  :  y  <  0.5  , 

(2.2.6b)  F (2)  :  a  >  0  ,  b  >  0 

and 

(2.2.6c)  F(3)  :  c  >  0  . 

Let  e  (Q)  =  r^  (Q)  -  exp  (Q)  be  the  error  for 
method  F  ^ 1 ^  . 

Definition 

If  the  free  parameters  are  chosen  such  that  the 
/  •  \ 

formulas  F^  '  are  made  exact  in  a  discrete  sense  (i.e., 
(Qq)  =  0  for  a  particular  Qq  =  hX^)  ,  then  the  for¬ 
mulas  F^  are  said  to  be  exponentially  fitted  at  Q  =  Qq  . 

By  definition,  we  require  the  following  relations 
to  hold  for  particular  Qq 1 s : 

(2.2.7a)  e (1) (Qq )  =  0  , 

e(2)  (Q0)  =  e(2)  (Qx)  =  0 


(2.2.7b) 
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and 


(2.2.7c) 

£  (3>  (Q„)  =  0 

where  QQ  and  are  obtained  from  two  different  eigen- 


values  Aq 

and  A^. 

The 

solution  of  each  equation  in  (2.2.7)  gives  the 

following 

values  for  the  parameters: 

(2.2 . 8a) 

y  =  -l/Q  -  1/  (exp (-Q)  -  1)  , 

(2.2.8b) 

a  =  2  (v(Q1)  -  v(Qq)  )/(QlV(Q0)  -  QQv  (Q^  )  , 

(2.2.8c) 

b  =  2(Qq-  Q1)/(QQv(Q1)  -  Q1v(Qq)) 

and 


(2.2. 8d) 

_  _  (Q2  +  6Q  +  12  -  exp  (Q)  (Q2  -  6Q  +  12 )J 

3  (exp (Q)  (Q2  -  2Q)  +  Q2  +  2Q) 

where 

v(z)=  (z2  (1  -  exp  (-z)  )  /  (  -  (2+z)  +  (2+z)  exp  (-z )  )  . 

We  know  that  conventional  methods  restrict  the  step- 
size  in  the  early  time  period.  The  main  advantage  of  this 
method  is  that  the  same  stepsize  used  in  the  later  time 
period,  is  used  during  the  early  time  period  of  integra¬ 
tion.  Another  advantage  is  that  the  free  parameters  can 
be  put  to  zero  after  the  transient  terms  do  not  contribute, 
thus  getting  higher  order  methods. 


n 
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2 . 3  Dahlquist 's  Smoothing  Technique 

Consider  the  differential  equation 

(2.3.1)  y  =  Ay  ,  Re(A)  <  0  ,  y(0)  =  1 

with  the  true  solution 


(2.3.2)  y(x)  =  exp(Xx)  . 

The  approximate  solution  of  (2.3.1)  at  xn=  nh 
obtained  using  the  trapezoidal  rule  is  yn=  (  (1+z) / (1-z) ) n 
where  z  =  hX/2 .  Clearly,  for  -z  very  large,  the  solution 
points  oscillate  about  zero.  To  damp  these  oscillations, 
Dahlquist  has  introduced  the  idea  of  smoothing.  Smooth¬ 
ing  consists  of  replacing  solution  points  yn  with  weighted 
averages  of  solution  points  in  the  neighbourhood  of  y  . 

For  example,  one  particular  smoothing  formula  is 

(2.3.3)  y  =  (y  ,  +  2y  +  y  ,-,)/4  . 

Jn  Jn-1  -*n  Jn+1  ' 


Dahlquist  suggests  using  this  formula  one  or  more  times 

in  calculating  the  solution.  To  define  a  typical  smooth- 

m 

ing  algorithm,  let  {n^}j_Q  be  any  sequence  of  integers 

satisfying  nn=  0,  n.  >  n.  ,  +  1  and  n  <  N,  where  m  is 
1  ^  0  j  j-1  m 

the  maximum  number  of  times  (2.3.3)  is  used  and  yn  /  yn  / 

.  .  .  ,  yn  represents  the  solution  points  at  which  smooth- 
m 

ing  is  desired.  The  smoothing  algorithm  for  (2.3.1)  is 
given  below.  The  input  consists  of  m,  n^,  n^, . . . ,n^,  N, 
z  and  y q • 


■ 
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1.  [READ]  READ  m,n^,  .  .  .  ,11^,  N,  z,  y  . 


2. 

[INITIALIZE] 

j 

1  ; 

y 

*<>• 

3. 

[TRAPEZOIDAL  RULE] 

y 

n  .  -1 

y (  (1+z)/  (1-z) )  3  . 

4. 

[SMOOTHING] 

y  «- 

(1  + 

2<l-z>  +  (i-z> 

5. 

[INCREMENT] 

j 

j  + 

1  ; 

IF 

j  * 

m  THEN  GOTO  3. 

6  . 

[FINISH] 

K  «- 

N  - 

n  +  1 . 
m 

7. 

[DECREMENT] 

K 

K  - 

1  ? 

IF 

K  = 

0  THEN  STOP . 

8. 

[UPDATE] 

y 

y (1+z)  /  (1-z) ; 

GO 

TO 

7. 

For  the  preceding  algorithm,  one  can  show  (Lindberg 
[17])  that  the  computed  solution  points  y^  satisfy 

(2.3.4)  yn  =  (  (l+z)/(l-z)  )n/(l-z2)m  =  yn/(l-z2)ra. 

It  is  important  to  notice  that  the  amplitude  of  the 

oscillations  of  y  is  smaller  than  those  of  y  ,  by  the 

u  n  u  n  -1 

2  m 

factor  (1-z  ;  .  However,  smoothing  methods  introduce 

new  errors.  For  the  above  algorithm,  the  error  e ^  is 
given  by 
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(2.3.5a)  ^(n)  =  exp(nhX)  -  yn/(l-z2)m  . 

The  error  of  the  trapezoidal  rule  defined  by 
(2.3.5b)  ^  (n)  =  exp(nhX)  -  y^ 

and  the  error  satisfy 

(2.3.6)  e1(n)  =  ^ (n)  +  yR(l  -  l/(l-z2)m)  . 

From  (2.3.6),  we  see  that  for  |z|<  1/  the  error 
is  smaller  than  .  However,  for  |z|  >1,  equation  (2.3.6) 
clearly  shows  that  the  amplitude  of  the  error  will  be 
growing.  We  note  that  as  m  increases,  the  amplitude  of 
the  oscillations  gets  smaller.  However,  Lindberg  shows 
that  for  the  non-stiff  components  (components  associated 
with  small  eigenvalues),  accuracy  is  lost  as  m  increases. 
Since  smoothing  itself  creates  new  errors,  this  error  is 
proportional  to  the  error  in  the  trapezoidal  rule  for 
linear  systems  with  constant  coefficients.  For  non¬ 
linear  systems  with  large  initial  values,  the  errors  may 
get  amplified. 

Lindberg  combines  the  smoothing  technique  with 
extrapolation  and  shows  that  without  extrapolation, 
smoothing  cannot  be  done  with  the  trapezoidal  rule.  We 
note  that  many  other  smoothing  formulas  exist,  but  only 
those  which  introduce  errors  commensurate  with  the  errors 
made  by  the  integration  formula  are  acceptable. 


■ 
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CHAPTER  III 


MULTISTEP  METHODS 


3 . 1  Gear's  Method 

For  the  system  (1.1.1),  consider  the  multistep 

method 


k 


k 

,  i  •  —  h  y  8 .  f  , 

n+l~D  ji0  :  n+1_3 


where  0  and  | |  +  |  3k  |  f  0 . 

For  stability  analysis,  we  are  interested  in  the 
behaviour  of  the  roots  of  the  polynomial 


(3.1.2) 


P  (O  +  hXj<j  (O  =  0 


9f/9y.  The  stability  domain  of  (3.1.1)  is  the  set  of 
hA  for  which  the  roots  of  (3.1.2)  are  inside  the  unit 
circle.  The  boundary  of  the  stability  domain  is  given 
by  the  locus  of  -p  (£)/a(£)  for  E,  =  exp(i0),  0  £  [0,2tt]. 

Gear  [11]  observed  that  as  h\  changes  from  0  to  -«> 
the  roots  of  (3.1.2)  move  from  the  roots  of  P(£)  to  thos 
of  a(£).  The  most  stable  a(£)  polynomial  of  degree  k 
is  cr(£)  =  £  •  By  using  the  stability  domains.  Gear  has 
shown  that  for  k  <  6,  these  methods  are  stiffly  stable. 
For  example,  if  k  =  6  then  6-^  satisfies  6^  <  -6.1  and  y 
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satisfies  y<  0.5.  Jain  and  Srivastva  [12]  have  shown  the 
existence  of  stiffly  stable  methods  upto  order  11. 

Using  the  corrector  formula  (3.1.1),  we  obtain  the 
fixed  point  iteration 


(3.1.3) 


(m+l)  _  _ ^  f  (x  ,  ,  y^jJ)  -  C 

n+1  M0  n+1  ^n+1 


where 

k 

C  =  y  [a-  y  , -,  .  +  h$.f  .]  . 

j=l  3  n+1-:  p:  n+1-3 

If  the  predicted  value  of  yR+^  is  y^| ,  this  fixed  point 
iteration  may  not  converge  for  large  values  of  h3f/3y, 
because  the  normal  convergence  condition  is  | h 3 ^ U f / Sy |  <  1. 

This  restriction  can  be  avoided  by  reformulating  the 
iteration  (3.1.3)  as 


(3.1.4a) 


(m)  (m+l) 

n+1 


t . y_r:v'  =  -hB0[f<y‘;i>-D-'yrb  -  c 


(m)  ,r  (m) 
n+1' 


where 


T (m)  =  [I  +  h3QD (ra) ] 


and  D  ^  is  a  suitably  chosen  matrix. 

(m) 


To  describe  how  D 
(m) 


should  be  chosen,  we  show  that 
(m) 


if  D  =  3f/3y  ,  evaluated  at  y^+^,  then  (3.1.4a)  is  a 
Newton  iteration. 

Dropping  the  subscripts,  we  can  write  (3.1.3)  as 


y  =  -h3Qf (y)  -  C  . 


. 
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Using  Newton's  method  to  find  y,  we  obtain 

(3.1.4b)  [I  +  h3Q9f (m)/3y]  (y  (m+1)  -  y (m) ) 

=  -(y  <m)  +  hgQf  (y  <m)  )  +  C)  . 

By  adding  to  both  sides  of  (3.1.4b)  and 

rearranging  terms,  we  obtain 

(3.1.4c)  T  (n>)  y  (m+1)  +  hgo  {y  (“+D  (3f  (“»/3y  -  D<m>) 

+  y  (m)  (D(m)  -  3f  (m)/3y) ) 

=  -h6o(f(y(m))  -  D(m)y(m))  -  C  . 

/  \ 

Choosing  the  matrix  D  ^  '  to  be  the  Jacobian  matrix  9f/3y, 

(m  \ 

evaluated  at  v  ,  yields  the  required  result. 

For  computational  purposes,  Gear  reformulates  equa¬ 
tion  (3.1.4a)  for  m  >  1  as  follows: 


(3.1.5)  y(m+1)  =  y(ra)-  P0[T(m,]"1(hf(y(m,)-d(m)) 

where 

d(m)  =  hf(y<ra-1))  +  hD(m-l)  (y(m)  _  ydn-Dj  _ 

Furthermore,  he  finds  that  given  d^t  the  sequence  {d 
satisfies 


(3.1.6)  d  (m+!)  =  d(m)  +  1(hf(y(m))  -  d (m)  )  . 


In  fact,  Gear  shows  that  (3.1.5)  and  (3.1.6)  can  be  used 
for  m  >  0  provided  that 


[  J 
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(0) 


(0) 


(y  +  ^  'ajyn+i-j  +  hsjyn+1_j))/e0 


Gear  [13]  has  shown  that  any  general  multistep  method 
(3.1.1)  can  be  written  as  a  multivalue  method.  Thus  if 
(3.1.5)  is  written  in  the  multivalue  form,  we  obtain 


(3.1.7a) 


(m+1) 


—n+1 


z  f  ( z  ) 

-n+1  -  —n+1  ' 


(3.1.7b) 

and 

(3.1.7c) 


z(0) 

-n+1 


=  Bz 


z  i-i 

—n+1 


-n 


=  z(M) 
-n+1 


where 


(m)  _  r  (m) 

-n  “  [yn  '  yn-l'  yn-2 ' *  *  * 'yn-k 


dnm'  '  hf n-1 ' *  *  * 'hfn-k] 


b  =  [-3  ,  0  ,  . . . ,  0  ,  1  ,  0  ,  . . .  0]  , 


F(2(rn),  =  IT(»)j-l(hf(yto))  -  d<*>) 


k-n+l 

and  the  matrix  B  is  given  by 


B 


. 
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The  coefficients  in  B  are  given  by 


nj  =  ~(aj  +  aj)//3o  ' 


h  -  -(6j +  Bj>/e0  - 


*  * 

where  ou  ,  3j  are  the  coefficients  associated  with  the 
following  predictor  formula 


r<0) 

n+1 


■*'*  *  "k 

1  (a .  y  ,  ,  .  +  h3 .  f  ,  . ) 

j=l  3  n+1-]  j  n+1-: 


Gear  reformulates  the  multivalue  method  (3.1.7), 
using  standard  techniques  (Gear  [12] )  into  a  multivalue 
formulation  using  vectors  an  of  scaled  derivatives  given 
by 


,h2y  /2 I 
i  n 


,  k-1  (k-1 )  /n  ,  .  ,  .t 
,h  y^  V  (k-1)  !  ] 


In  the  Nordsieck  formulation,  the  multivalue  method  has 
the  form 


(3.1.8a) 


a  |  =  A  a 
—n+1  — n 


(3.1.8b) 


(m+i)  _  <m)  G(a<™>) 

—n+1  —n+1  —  —n+1 


and 

(3.1.8c) 


a  =  a(M) 
-n+1  -n+1  ' 


' 


’it*  -  '  . 
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Initially,  all  the  unknown  derivatives  of  a^  are 
set  to  zero.  Then  given  a  ,  Gear's  method  uses  equation 
(3.1.8a)  to  predict  a^|.  The  corrector  equation  (3.1.8b) 
is  then  iterated  without  changing  the  Jacobian  matrix 
until  the  number  of  iterations  M  equals  3.  If  convergence 
is  not  obtained,  the  Jacobian  matrix  8f/9y  is  evaluated 
again.  If  this  happens  a  second  time,  the  step  is  reduced 
by  one-half  and  the  process  starts  again.  An  estimate  of 
the  local  error  is  obtained  by  taking  the  difference  of  the 
predicted  and  corrected  values  at  each  step.  The  stepsize 
is  increased  or  decreased  by  a  factor  of  2  depending  on  the 
error  tolerence.  Fcr  changing  the  order  k,  the  maximum 
step  for  orders  k-1,  k  and  k+1  is  calculated.  Whenever 
the  stepsize  or  the  order  is  changed,  k  steps  must  be  taken 
before  we  can  change  either  the  stepsize  or  the  order. 


. 

. 

' 
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3 . 2  Fourth  Order  Exponentially  Fitted  Method 

Liniger  and  Willoughby  [19]  introduced  the  concept 
of  exponential  fitting.  Recently  Bjurel  et  al.  [2]  have 
defined  multistep  methods  which  are  exponentially  fitted 
in  the  complex  hX-plane.  Liniger  and  Willoughby  [19] 
have  developed  exponentially  fitted  methods  of  order  three. 
In  this  section,  we  develop  a  new  fourth  order  exponen¬ 
tially  fitted  method. 

(4) 

Consider  the  following  integration  formula  F  of 
order  p  =  4 


(3.2.1a)  F(4>:  Yn+2  =  Yn+1  +  <A*n  +  B*n+l+Cyn+2) 


+  240  +  Eyn+l+FW  +  T 


where 


A  =  50a  +  26b  -  65  , 


B  =  -40a  -  52b  +  220  , 


C  =  -10a  +  26b  +  85  , 

(3.2.1b) 

D  =  20a  +  13b  -  30  , 


E  =  40a  , 

F  = -13b 


and  a ,  b  are  real  parameters.  The  truncation  error  T  is 
given  by 


*  1$ 

■ 

s. 
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(23  -  13b-  10a)  ,5  (5)  ,,  6, 

T - 1440 -  h  Yn  +  0  (h  > 


showing  that  (3.2.1a)  is  a  fourth  order  method  at  least. 

(4 ) 

For  a=b=l,  formula  F  reduces  to  a  sixth  order  strongly 
stable  formula  given  by 

Yn+2  =  Yn+1 +  2T0  (llYn  +  128Yn+l  + 

+  210  (3Yn  +  4°Yn+l  '  13Yn+2>  * 

For  a  =  1/2  and  b  =  20/13,  formula  F  ^  reduces  to  the 
well  known  formula  associated  with  the  second  diagonal 
Pade  approximation  given  by 


Y 


n+1 


Yn  +  I  <Yn+l+  V  "  II  <Yn+l  '  V  ' 


(4  ) 

The  multistep  method  F  ,  when  applied  to  the 
scalar  equation  y  =  Ay,  Re (A)  <  0,  yields  the  following 
difference  equation 


(3.2.2a)  yn+2  (2  40  -  CQ  -  FQ2)  +  yn+1  (-2  40  -  BQ  -  EQ2) 


+  yn  (-AQ  -  DQ  )  =0 


where  Q  =  hA.  Its  characteristic  polynomial  x(z'"Q)  is 


X (z;Q)  =  z2  (240  -  CQ-  FQ2)  +  z (-240  -  BQ  -  EQ2 ) 


-  AQ  -  DQ' 


(3.2.2b) 


' 
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The  polynomial  has  two  roots,  both  of  which  must 
be  less  than  1  in  modulus.  The  principal  root  should 
approximate  exp (Q) .  The  most  stable  choice  for  the 
parasatic  root  is  zero.  Thus,  we  choose  the  parameters 
in  x ( z ; Q)  so  that 


(3.2.3a)  x(exP(Q);Q)  =  0 

and 

(3.2.3b)  x(0/'  Q)  =  0  . 


Solving  (3.2.3)  for  a  and  b,  we  obtain  for  Q  ^  -2.5 


and 


where 


and 


a  =  (65  +  30Q  -13(2+  Q)S(Q))/(50  +  20Q) 

b  =  S(Q) 

S(Q)  =  N(Q)/D(Q)  , 

N(Q)  =  (5  +  2Q)  (exp(Q)  (85Q  -  240)  +  220Q+  240) 

-  Q (65  +  30Q) p (Q)  , 

D  (Q)  =  13Q((5  +  2Q)  (exp  (Q)  (Q  -  2)  +  4)  -  (2  +  Q)p(Q)) 

p  (Q)  =  exp(Q)  +  4(1  -  Q)  . 


Since  a  and  b  are  chosen  to  make  the  parasatic  root 


zero,  the  principal,  root  of  (3.2.2b)  is 


. 
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240  +  BQ  +  EQ2 

240  -  CQ  -  FQ2 


(4 ) 

To  prove  that  F  is  A-stable,  we  use  the  follow¬ 


ing  theorem  due  to  Liniger  and  Willoughby  [20] : 

Theorem 

Let  r  (Q)  be  a  non-constant  analytic  function  for  Real 
Q  <  0.  Then  |r(Q) |  <  1  for  Real  Q  <  0  if  and  only  if  both 
of  the  following  are  true: 

(3.2.4a)  jr(Q) |  <  1  on  Real  Q  =  0 


and 


(3.2.4b)  Lim  Ir  (Q)  |  <  1  . 

I  Q  I  -*  00 


Theorem 


The  necessary  and  sufficient  conditions  for  the  for- 
(4 ) 

mula  F  to  be  A-stable  are 


(3.2.5a)  a/b  <  13/40 


(3.2.5b)  I C | 


< 


(3.2.5c) 


o  <;  c 


and 


( 3 . 2 . 5d) 


0  <  b 


where  B  and  C  are  defined  by  (3.2.1b). 
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Proof 

(4 ) 

F  is  A-stable  if  and  only  if  |r(Q) |  <  1  for 
Real  Q  <  0.  First  we  show  that  if  |r(Q)  j  <1,  then  con¬ 
ditions  (3.2.5)  hold. 

Let  Q  +  on  the  real  axis.  Then  we  must  have 

|E |  <  j  F  j  ,  which  using  (3.2.1b)  yields  equation  (3.2.5a). 

Since  |r(Q) |  <  1  for  Real  Q  <  0,  all  the  poles  of 

r  (Q)  are  in  the  right  half  complex  plane.  Thus  the  zeros 

2 

Zj  of  240  -CQ-FQ  ,  satisfy  Real  z^  >  0,  which  is  true 
if  and  only  if  F  <  0  and  C  >  0.  Hence  (3.2.5c)  and 
( 3 . 2 . 5d)  hold. 

Finally ,  | r (Q)  |  <  1  implies  | B  +  EQ  j  <  | C  +  FQ |  . 

Letting  Q  0  gives  conditions  (3.2.5b). 

To  show  that  these  conditions  imply  that  | r (Q) |  <  1 
for  Real  Q  <  0,  note  that  conditions  (3.2.5c)  and  (3.2.5d) 
imply  that  r (Q)  is  analytic  in  the  left  complex  plane. 

We  finally  note  that  conditions  (3.2.5a)  and  (3.2.5b) 
insure  that  the  conditions  of  the  previous  theorem  are 
satisfied. 


Q .  E  .  D  . 


* 
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CHAPTER  IV 

TEST  PROBLEMS  AND  CONCLUSIONS 

4 . 1  Introduction 

In  this  chapter,  results  of  test  problems,  using 
both  the  Liniger  and  Willoughby  method  and  the  new  fourth 
order  method  are  given.  The  test  problems  chosen  repre¬ 
sent  both  linear  and  non-linear  systems,  with  moderate 
and  highly  stiff  combinations.  Since  | h 3 f / 3 y |  is  large 
for  stiff  systems,  Newton's  method  is  used  for  solving 
the  implicit  difference  equations,  which  arise  in  apply¬ 
ing  each  method.  The  inputs  to  each  computer  program 
were  the  initial  values,  eigenvalues  at  x  =  0,  the  error 
tolerance  and  the  Jacobian  8f/3y.  For  systems  whose 
Jacobian  is  a  function  of  time,  we  used  the  eigenvalues 

found  for  x  =  0.  The  first  solution  point  for  the  fourth 

(3 ) 

order  method  was  computed  by  formula  F  ,  using  the  step- 
size  h/8,  where  h  is  the  stepsize  for  the  fourth  order 
method . 


. 


' 
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4 . 2  Computational  Equations  for  the  Liniger  and 
Willoughby  Method 

Applying  formula  to  (2.1.1),  the  following 

difference  equation  is  obtained 

(4-2-1)  Yn+1  -  (1-v)hf„+l  “  Yn  +  Hhfn  • 


Equation  (4.2.1)  is  solved  by  the  following  Newton 


iteration 

(4.2.2) 

(P)  _ 
n+1 

where 

A (P>  _  V(P+1)  _ 

v(p) 

Vl  "  Yn+1 

Yn 

JnA  m  3f(xn+l' 

y(p> 

yn+l 

f (P )  _  f (x 

n+1  1  n+1' 

y(p|) 

Yn+1 

with  = 

y  and  y  =  y 

7  n  Jn  2 

(M) 
n  * 

the  maximum 

number  of  iterations 

where 

y(P+D  .  (P) 

yn+l  yn+l 

<  e 

n 


n 


n+1  -*n+l 


with  e  a  fixed  tolerence. 

(2) 

Applying  formula  Fv  '  to  (2.1.1),  the  following 
difference  equation  is  obtained 


40 


yn+i  -  h  (1+a)  f n+1/2  +  h2  (b+a)  Jn+lW4 

-  y  -  h (1-a) f  /2  -  h2  (b-a) J  f  /4  =  0  . 

J  n  iv  n  iv 

By  using  the  Newton  iteration,  we  obtain 
(4.2.3)  (41  -2h(l+a)J^P|  +  h2  (b+a)  [ ]  2 ) 

=  4(^n-  yn+i)+h2(2(1+a)I-h(b+a)Jn?i)fn+i 

+  2h ( 1-a) f  +  h2  (b-a) J  f 
n  n  n 

( 3 ) 

The  computational  equation  for  F  follows  from  (4.2.3) 
by  substituting  b  =  1/3  and  c  =  a. 


■ 


' 
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4 . 3  The  Computational  Equation  for  the  Fourth  Order 

Exponentially  Fitted  Method 


Applying  formula  (3.2.1a)  to  (2.1.1),  the  follow¬ 


ing  difference  equation  is  obtained 


+  h2 (DJ  f 
n  n 


+  EJ  +  FJ  ^0f  „ )  /240 

n+1  n+1  n+2  n+2 ' 


where  the  coefficients  A,  B,  C,  D,  E  and  F  are  defined  by 
(3.2.1b) . 

Following  Liniger  and  Willoughby's  technique  of 
linearizing  the  approximations  (see  section  4.2),  the 
computational  equation  becomes 


(2401  -  hCJ^  -  h2F  [J 


+  h  (AI  +  hDJ  )f  +  h(BI  +  hEJ 


n  n 


n+1  ^  ^  n+1 


+  h (Cl  +  hFJ 


(P)  x f (P) 
n+2'  n+2 


where  y^°i  =  y  , ,  is  used  for  starting  each  iteration. 
n+2  n+l 


n+2 
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4 . 4  Ehie's  Linear  Problem  I  [8] 

(4.4.1)  y  =  -  y  +  95z  ,  y(0)  =  1  , 

z  =  -  y  -  9 5z  ,  z  (0 )  =  1  ,  x  e  [0,1]. 

Eigenvalues  of  the  Jacobian  at  x  =  0  :  -96  ,  X^=  -2. 

True  solution:  y  =  (95  exp(-2x)  -  48  exp (-96x) ) /47 

z  =  (48  exp(-96x)  -  exp(-2x))/47  . 


Table  2 


The  Solutions  of  Equation  (4.4.1) 


STEP 

METHOD 

yd) 

( | ERROR | ) 

z (1)  x 102 

( | error] xio2) 

0.0625 

f(3) 

0.2735430 

(7xl0“6) 

-0.2879400 

( 7xl0~6 ) 

f(4) 

0.2735503 

(3xl0-7) 

-0.2879477 

(4xl0_7) 

0.03125 

f(3) 

0.27354952 

(5xl0~7 ) 

-0.28794687 

(5xl0_7) 

f(4) 

0.27355005 

( lx 10~  8 ) 

-0.28794742 

(lxlO-8) 

0 . 015625 

f(3) 

0.27355000 

( 4x 10_  8 ) 

-0.28794737 

( lxl0~  8 ) 

TRUE  SOLUTION 

0.27355004 

-0 .28794741xl0“2 

( " 
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4 . 5  Ehle's  Linear  Problem  II  [8] 

(4.5.1)  y  =  -  y  +  95z  ,  y(0)  =  1  , 

z  =  -  y  -  95z  ,  z  (0 )  =  -1/95  ,  xe  [0,1]. 

Eigenvalues  of  the  Jacobian  at  x  =  0  :  X^=  -96,  X2=  -2. 

True  solution:  y  =  exp (-2x) 

z  =  -  (exp  (-2x) ) /95  . 


Table  3 


The  Solutions  of  Equation  (4.5.1) 


STEP 

METHOD 

yd) 

( I  ERROR  I  ) 

z (1)  x  102 
( | ERROR | xlO2) 

0.0625 

p(3> 

0.135331816 

(3xl0“6) 

-0.14245454 

(3xl0-6) 

f(4) 

0.1353354305 

(lxlO-7) 

-0.1424583479 

(lxlO-7) 

0.03125 

f(3) 

0 . 1353350304 

( 2x 10~  7 ) 

-0.14245792 

(2xl0~7) 

f(4) 

0 . 1353352887 

(5xl0“9) 

-0.1424581987 

(5xl0_9) 

0.015625 

f(3) 

0.1353352666 

( lx 10_  8 ) 

-0.14245817 

( 2x 10~  8 ) 

TRUE  SOLUTION 

0.1353352832 

-0.1424581928x1 0~2 

- 


- 
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4 . 6  Fowler  and  Warten's  Linear  Problem  [9] 

(4.6.1)  y  =  -2000y  +  lOOOz  +  1  ,  y(0)  =  0  , 

z  =  y  -  z  ,  z(0)=0f  x  e  [0,4]  . 

Eigenvalues  of  the  Jacobian  at  x  =  0  :  X^=  -2000.5,  A2=-0.5. 


Table  4 


The  Solutions  of  Equation  (4.6.1) 


STEP 

METHOD 

y (4)  X 103 
( | ERROR | xlO  3 ) 

z (4)  x  103 
( | ERROR | xlO3) 

0.0625 

f(3) 

0.93226472 

(6xl0~8) 

0 . 864563300 

(lxlO"7} 

f(4) 

0.93226466 

(EXACT) 

0 . 8645631887 

(lxlO-9) 

0.03125 

f<3> 

0.93226467 

(lxlO-8) 

0 . 864563203 

(lxlO-8) 

f(4) 

0.93226466 

(EXACT) 

0 . 8645631898 

(lxl0_1°) 

0.015625 

f<3> 

0.93226466 

(EXACT) 

0 . 864563191 

(2xl0-9) 

TRUE  SOLUTION 

R-K  4  RALSTON 
h  =  0.001 

0.93226466xl0“3 

0 . 8645631899xl0-3 

. 


. 


' 

, 

4.7  Variable  Eigenvalue  Linear  Problem  [10J 


(4.7.1)  Y  =  AY  ,  Y ( 0 )  =  [2/0]t  ,  x  e  [0,1] 

where 


1+c 

1-c 

1-c 

1+c 

Eigenvalues  of  the  Jacobian  at  x  =  0  :  X^=  -c  ,  X2=  -1 

True  solution:  y1  =  exp(-x)  +  exp(-cx) 

y2  =  exp(-x)  -  exp(-cx)  . 

We  consider  the  particular  case  c  =  16 . 


Table  5 


The  Solutions  of  Equation  (4.7.1) 


STEP 

METHOD 

Yp  U) 

( | ERROR | ) 

y2  (1> 

( | ERROR | ) 

0.0625 

f(3) 

0.367879437 

(IxlO-7) 

0 . 36787921266 

( lxl0~  7 ) 

f(4) 

0.3678795562 

(  2xl0"9) 

0.36787933114 

(2xl0“9) 

0.03125 

f(3> 

0.36787954641 

(7xl0_9) 

0.36787932134 

(7xl0-9) 

f(4) 

0 . 36787955378 

(8xl0_11) 

0.3678793287 

(lxl0_1°) 

TRUE  SOLUTION 

0.36787955370 

0.3678793286 

••  fr 

.r  ;  > 


% 


V 


C  J  sr  - 

*NNM to*,  j  •  ..  **  . 


*r 

-  ««i  ..*.  -  - 


*V*«44 
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4 . 8  Lapidus 1 s  Linear  Problem  [15] 

(4.8.1)  Y  =  AY  ,  Y ( 0 )  =  [2/l/2]t  ,  x  e  [0,1] 


where 


A 


-0.1  -49.9  0 

0  -50  0 

0  70  -120 


Eigenvalue  i  of  the  Jacobian  at  x =  0  :  Xn=  -120, 


X3=  -0.1. 

True  soliii  .on: 

^2 


y 


3 


exp  (-0. lx)  +  exp(-50x) 
exp  (-50x) 

exp(-50x)  +  exp(-120x)  . 


-50, 


Table  6 


The  Solutions  of  Equation  (4.8.1) 


STEP 


0.0621 


0.031: 


METHOD 

Yld) 

( | ERROR | ) 

yk(i)  k=2^3 

(  | ERROR |  ) 

f(3) 

0.904837417863 

(2xl0_1°) 

0.14xl0"26 

(0.19xl0-21) 

p(4) 

0.90483741803630 

(4xl0-13) 

0 .25xl0-17 

(0.25xl0“17) 

f<3> 

0.904837418022840 

(lxlO-11) 

O.llxlO-21 

(0.8xl0_22) 

f(4) 

0.9048374180359725 

-14 
(1x10  xq) 

0 . 31xl0-21 

(O.lxlO-22) 

SOLUTION 

0.9048374180359595 

0.19xl0"21 

- 


■ 

47 


4 . 9  Fowler  and  Warten 1 s  Non-linear  Problem  [9] 

(4.9.1)  y  =  -lOz  +  3000  (1-y2)  ,  y(0)=  1, 

z  =0.04  (1-z)  -  (1-y)  z  +  10~4  (1-y)  2  ,  z(0)=  1,  xe  [0,2.6] 

Eigenvalues  of  the  Jacobian  at  x=  0  :  X^=  -6000,  X2=  -0.04. 


Table  7 


The  Solutions  of  Equation  (4.9.1) 


STEP 

METHOD 

y  (2.6) 

1  (ERROR)  I 

z  ( 2 . 6  )  x  10 
( | ERROR | x 10) 

0.1625 

f(3) 

0.999939925473 

(EXACT) 

0.9876685530 

(4xl0~9) 

f(4) 

0.999939925473 

(EXACT) 

0.98766854924 

(lxlO-11) 

0.08125 

f(3) 

0.999939925473 

(EXACT) 

0.9876685497 

(5xl0_1°) 

f(4) 

0.999939925473 

(EXACT) 

0.98766854924 

(lxlO-11) 

TRUE  SOLUTION 

R-K  4  RALSTON 
h  =  0.00040625 

0.999939925473 

0.98766854923x1 0_1 
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4.10  Lawson's  Non-linear  Problem  [16] 

(4.10.1)  y  =  (-i+z2)y+  (l+z)z  ,  y(0)=-l 

z  =  -y  +  (-19+2y+z^)  z  ,  z  (0  )  =  1  ,  x  e  [0 , 1]  . 

Eigenvalues  of  the  Jacobian  at  x=  0  :  -19.9,  X2=  -0.1. 


Table  8 


The 

Solutions 

of  Equation  (4.10.1) 

yd) 

z  (1)  x  10 

STEP 

METHOD 

( | ERROR | ) 

( | ERROR | x 10) 

f(3) 

-0 . 33062729 

0.17849359 

0.03125 

(4xl0-6) 

(2xl0_6) 

f(4) 

-0.33063158 

0.17849585 

(7x10  ') 

(3x10  ') 

f(3) 

-0 . 33063060 

0.178495338 

-7, 

-7, 

0.015625 

(3x10  ') 

(2x10  ') 

f(4) 

-0.33063088 

0.17849548 

(3xl0-8) 

(lxlO-7) 

TRUE  SOLUTION 

-1 

R-K  5 

-0.33063085 

0.1784955x10 

LAWSON  [16] 
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4.11  Gear  1 s  Non-linear  Problem  [11] 


(4.11.1) 

y  =  -(55  + 

z)y  +  65z  , 

y (0)  =  l  , 

z  =  0.0785 

(y-z)  r 

z(0)  =  1  , 

u  =  0 . 1  y  , 

u (0)  =  0  ,  x  £  [0 ,1] 

Eigenvalues  of  the  Jacobian  at  x=  0  :  X^=  -55, 
X2  ,X3  =  0.0625  ±  O.Oli  . 


Table  9 


The  Solutions  of  Equation  (4.11.1) 


STEP 

METHOD 

yd) 

( | ERROR | ) 

z(l) 

( | ERROR | ) 

U(l) 

( I  ERROR  I ) 

0 . 125 

f(3> 

1.1955191535 

(8xl0“9) 

1.013991515 

(7xl0-9) 

0.118510694 

(lxl0“8) 

f<4> 

1.19551914432 

(9xl0_1°) 

1.01399150753 

(8xl0_1°) 

0.1185106830 

(lxl0~9) 

0.0625 

f(3> 

1.19551914618 

(lxl0“9) 

1.0139915091 

(8xl0_1°) 

0.1185106852 

(lxl0“9) 

f(4> 

1.19551914499 

1.01399150814 

0.1185106838 

(2xl0_1°) 

(2xl0_1°) 

(2xl0_1°) 

TRUE  SOLUTION 
GEAR ' S 

PROGRAM  [12] 

1.195519145189 

1.01399150830 

0.1185106840 

. 
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4 . 12  Liniger  and  Willoughby's  Non-linear  Problem  [19] 

(4.12.1)  y  =  a1a5~  (a1+a2  (y+a4)  (y+1)  )  p  (y  ,  z)  ,  y(0)=0 

z  =  -  (1  +  a^z2  )  p  (y ,  z )  ,  z(0)  =  0,  xe  [0,100]. 

where 

P  (y  ,z)  =  a5  +  y  +  z  , 

a-^=a,-=a2=  1  ,  a^=  1000  and  a^=  0.01  . 

Eigenvalues  of  the  Jacobian  at  x=  0  :  A  =  -1012,  -0.01. 


Table  10 


The  Solutions  of  Equation  (4.12.1) 


STEP 

METHOD 

y  (loo ) 

(Difference  between 
F4  (h-0 .06  25)  and  yn| 

z  (100) 

| ERROR | 

0.125 

f<3) 

-0.99163973 

(2xl0~6) 

0.9833336 

(3xl0~6 ) 

f(4) 

-0.991641966 

(8xl0-8) 

0.983336240 

(lxlO-7) 

0.0625 

f(3) 

-0.991640918 

(lxlO-6) 

0.983335044 

(lxlO-6) 

f<4> 

-0.991642044 

(ASSUMED  EXACT) 

0.983336329 

(3xl0“8) 

TRUE  SOLUTION 

R-K  4 

ODEN  [21] 

NOT  REPORTED 

0.983336361 

. 
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4.13  Moore's  System  [7] 

(4.13.1)  y  =  -lOOz  +  lOOv  -  c(z-v)2  ,  y(0)  =  1  , 

z  =  y  /  z (0)  =  o  , 

u  =  (100 z -  10 Iv  +  c (z-v) 2-  dv2)/10  0  ,  u(0)  =  0  , 
v=u,  v  (0)  =  0  ,  x  6  [0,65]  . 

where  c  =  d  =  0  for  the  linear  case. 

Eigenvalues  of  the  Jacobian  at  x  =  0  :  L  =  -0.5 ±  lOi 

J-  /  2 

A  ^  ^  =  -0.005  ±  0 . li . 

Table  11 

The  Solutions  of  Equations  (4.13.1) 


c  =  d  =  0 


STEP 

METHOD 

y  (65) xlO2 

z  (65) xlO1 

u  (65) xlO2 

v(65) xlO 

(  | ERROR | ) 

0.0625 

f(3) 

0.70 

0.12 

0.70 

0.126377 

(7xl0_1) 

f(4) 

0.70 

0 . 12 

0.70 

0 . 126377 

(7xl0_1) 

TRUE  SOLUTION 
R-K  4 

ODEN  [21] 

NOT 

REPORTED 

NOT 

REPORTED 

NOT 

REPORTED 

0.196403 

. 
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4.14  Conclusions 

We  have  discussed  those  methods,  which  we  regard  as 
being  the  most  important  for  solving  stiff  ordinary  dif¬ 
ferential  equations.  There  are  many  other  methods  in  the 
literature  [Bjurel  et  al.  [2]),  but  these  appear  to  be 
too  special  in  the  sense,  that  they  are  applicable  only 
to  problems  satisfying  very  special  constraints.  In  fact, 
even  those  methods  surveyed  in  this  thesis  are  not  appli¬ 
cable  in  general. 

Of  the  single  step  methods,  the  modified  R-K  methods 
simply  increase  the  stability  boundary.  Calahan's  method 
and  Lawson's  method  for  linear  systems  appear  to  be  most 
promising  among  R-K  type  methods.  Dahlquist's  smoothing 
technique  is  quite  effective  for  damping  the  oscillations 
in  linear  systems  with  constant  coefficients.  Experimental 
results  reported  by  Lapidus  et  al .  [15]  indicate  that  the 

Liniger  and  Willoughby  method  is  the  best  among  single  step 
methods,  although  Calahan's  method  is  computationally 
faster  for  moderately  stiff  systems.  Also,  it  was  found 
that  for  Moore's  system  (4.13),  which  has  large  eigen¬ 
values  near  the  imaginary  axis,  the  Liniger  and  Willoughby 
method's  performance  was  quite  poor. 

The  most  promising  multistep  methods  seem  to  be 
Gear's  method  and  Dahlquist's  SAPS  method  [7].  Oden  [21] 
reports  that  for  the  SAPS  method,  the  number  of  function 
evaluations  are  less  than  that  of  the  Gear's  method.  We 
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note  that  Oden  has  chosen  only  those  stiff  systems,  which 
satisfy  conditions  imposed  by  the  SAPS  method.  Further¬ 
more,  Gear's  method  is  applicable  to  both  stiff  and  non¬ 
stiff  systems  and  it  has  both  automatic  stepsize  and 
variable  order  control.  Thus  among  multistep  methods, 
we  would  prefer  to  use  Gear's  method. 

The  second  objective  of  this  thesis  was  to  experiment 
with  the  Liniger  and  Willoughby  methods.  In  the  process 
of  experimentation,  we  developed  a  new  fourth  order 

method.  Since  the  new  method  involves  extra  function 

( 3 ) 

evaluations,  the  cost  of  n  steps  for  both  the  F  and 
(4 )  . 

F  methods,  is  given  in  Table  12.  Here  M  represents 
the  maximum  number  of  iterations  per  step  and  £  represents 
the  dimension  of  the  system.  We  assume  that  the  overhead 
is  the  same  for  both  methods. 


Table  12 

The  Cost  of  n  Steps  for  F^  and  F^ 


TYPE  OF 
STIFF  SYSTEM 

METHOD 

MATRIX 

INVERSION 

MATRIX . 
SCALAR 

MATRIX. 

VECTOR 

MATRIX- 

ADDITION 

VECTOR- 

ADDITION 

NON-LINEAR 

y  =  f  (y) 

f(4> 

nM£3 

2 

9nM£ 

3nM£2 

9nM£2 

6nM£ 

f(3) 

rmz3 

2 

6nM£ 

2nM£2 

5nM£2 

6nM£ 

LINEAR 

WITH 

CONSTANT 

COEFFICIENTS 

f(4) 

£3 

£2 

n£2 

0 

8n£ 

f(3) 

£3 

£2 

n£2 

0 

5n£ 
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For  the  non-linear  systems,  we  see  that  the  cost 

(4 ) 

of  F  is  higher  in  the  columns  of  MATRIX . SCALAR  and 
MATRIX-ADDITION.  For  linear  systems  with  constant  co¬ 
efficients,  we  note  that  the  cost  of  n  steps  for  F  ^ 
is  more  in  the  VECTOR- ADDITION  column  only.  In  fact, 

from  Tables  2-6,  we  obtain  that  F^^'s  performance  with 

•  •  ( 3 ) 

a  stepsize  2h  is  the  same  as  the  performance  of  F 

with  a  stepsize  h.  To  check  the  tradeoff  of  high  cost 
( 4 ) 

of  F  for  non-linear  systems,  the  execution  times  of 
both  methods,  for  the  non-linear  highly  stiff  system 
(4.12)  are  given  in  Table  13.  Here  error  represents  the 
maximum  error  in  either  y(100)  or  z(100). 


Table  13 

The  Excess  Time  of  F^  Over  F^ 


STEP 

METHOD 

! ERROR | 

EXECUTION 
TIME  IN 
SECS. 

EXCESS 

TIME 

0.125 

f(3) 

3xl0-6 

28.86 

2.40  secs. 

f(4) 

lxio"7 

31.26 

0.0625 

f(3) 

lxlO-6 

64.36 

1.75  secs. 

f<4) 

3xl0“8 

66.11 

' 


' 
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Thus  we  see  that  even  with  h  =  0.0625,  F  ^ does 

( 4 ) 

not  yield  the  accuracy  as  is  attained  by  F  with 
h  =  0.125.  We  note  that  the  increase  in  execution  time 
for  h  =  0.0625  is  just  2.7%,  while  the  error  is  reduced 
by  two  digits  at  least. 

In  conclusion,  noting  the  high  cost  and  higher 
(4 )  (4 ) 

accuracy  of  F  ,  we  would  say  that  F  should  be  used 

for  linear  systems  with  constant  coefficients.  For  non- 

( 4 ) 

linear  systems,  F  should  only  be  used,  if  high  accuracy 
is  desired. 
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